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We present general arguments for the importance, or lack thereof, of the structure in the charge 
distribution of counterions for counterion-mediated interactions between bounding symmetrically 
charged surfaces. We show that on the mean field or weak coupling level, the charge quadrupole 
contributes the lowest order modification to the contact value theorem and thus to the intersur- 
face electrostatic interactions. The image effects are non-existent on the mean-field level even with 
multipoles. On the strong coupling level the quadrupoles and higher order multipoles contribute 
additional terms to the interaction free energy only in the presence of dielectric inhomogeneities. 
Without them, the monopole is the only multipole that contributes to the strong coupling electro- 
statics. We explore the consequences of these statements in all their generality. 



I. INTRODUCTION 

The assembly of colloidal building blocks with designer 
engineered size, shape, and chemical anisotropy seems to 
be the next step in the fundamental and applied colloid 
science and science of soft materials, vigorously pursued 
by many researchers [l], 0]. Moving beyond the tradi- 
tional systems implies the creation of designer colloidal 
building blocks such as colloidal molecules, janus spheres, 
and other patchy particle motifs that have very differ- 
ent propensities for self-assembly 0]. In order to real- 
ize this goal one needs detailed control over various as- 
pects of colloid geometry at the nanoscale {e.g., aspect ra- 
tio, faceting, branching, roughness) and microscalc {e.g., 
chemical ordering, shape gradients, unary and binary col- 
loidal "molecules"). The ability to spatially modify the 
surface structure of colloids with designed chemical het- 
erogeneity, e.g., to form a patchy surface structure, seems 
to be becoming a realistic goal. Fabrication of stable 
anisotropic microcapsules was recently accomplished by 
the layer-by-layer polyelectrolyte adsorption technique 
combined with particle lithography technique to pro- 
duce anisotropic polymer microcapsules with a single 
nanoscale patch Granick and co-workers Q recently 
reported a highly scalable synthetic pathway for creating 
bipolar janus spheres, i.e., particles that consist of op- 
positely charged hemispheres. Such colloidal "animals" 
which are probably the simplest example of patchy col- 
loids, exhibit orientation-dependent interactions that go 
together with localized patches of like/unlike charges or 
hydrophobic/hydrophilic regions. This heterogeneous in- 
teraction landscape promotes the formation of larger col- 
loidal molecules and clusters that are themselves patchy. 

These advances in the nano- and microscopic tailor- 
ing of (charged) colloids motivated various approaches 



to generalizations of the existing theories of electrostatic 
interactions in charged colloids by explicitly including 
the structure of the counterions as embodied by their 
multipolar moments. The inclusion of structured coun- 
terions of a dipolar 0, Q and quadrupolar 0, H| type 
into the theory of electrostatic colloidal interactions has 
brought fourth some of the salient features of the coun- 
terion structure effects, which on face value appear to 
be quite distinct from the standard Poisson-Boltzmann 
framework. In view of these advances in the study of 
interactions between charged colloids it thus seems ap- 
propriate to explore the ramifications of the emerging 
paradigm if applied to these more complicated structured 
colloidal molecules. Especially the non-spherically sym- 
metric charge distribution of microscopically tailored col- 
loidal particles might lead to some unexpected properties 
of electrostatic interactions in this type of systems. With 
this in mind we thus embark on a thorough examina- 
tion of the consequences of multipolar charge distribution 
of mobile counterions that mediate interactions between 
charged (planar) macroions. 



II. MODEL 

In what follows we will consider a system of fixed 
macroions of surface charge distribution po (r) in an aque- 
ous solution, described as a dielectric continuum with 
a dielectric constant s and temperature X, containing 
N neutralizing counterions. Counterions are assumed to 
be pointlike particles but they do posses a rigid internal 
structure described by a charge distribution p{r;ILi,uJi) 
that we assume can be written as a standard multipolar 
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FIG. 1: Geometry of the system. Two charged walls with dif- 
ferent dielectric constants and structured counterions in be- 
tween. Counterions are treated as pointlike particles but with 
internal uniaxial structure that can be viewed, for example, 
as uniformly charged rods. 



expansion 



9] 



p(r; R^, Wi) = e qS(r - R 4 ) - p (n; • V)8(r - Rj 
+t (ru- V) 2 5(r-R 4 ) + -- - , 



(1) 



for the i-th counterion located at R^ with u)i being the 
orientational variables specifying the angular dependence 
of the counterion charge density. The monopolar mo- 
ment of each counterion is eoq, where q is the charge 
valency and eo the elementary charge, p = po n is the 
dipolar moment and Q = to n ® n the quadrupolar 
moment with director unit vector n. We must empha- 
size here that the quadrupolar expansion in (J]) is not 
general but adequately describes only a uniaxial counte- 
rion, e.g., a charged particle of rod-like structure, Fig. [1] 
One possible implementation of counterions possesing 
just monopolar and quadrupolar moment is a uniformly 
charged rod with charge eoq and length I. Counterion's 
quadrupolar moment is then to = eoql /24. Another 
possibility is a negative charge (—e^) in the center and 
two positive charges (+ei) located at both ends of the 
rod leading to monopolar and quadrupolar moments of 
eoq = 2ei — 62 and to = e\l 2 , respectively. 

Counterions interact via a Coulomb interaction poten- 
tial u(r, r') so that the interaction energy of two given 
counterions i and j at set positions and orientational con- 
figurations is obtained by integrating the Coulomb inter- 
action over their internal orientational degrees of freedom 
as 



J J dr'drp(r; Rj, «i)it(r, r')p(r'; R^ , uj 



(2) 



Although our formalism in this section is in general ap- 
plicable to macroions of arbitrary shape and charge dis- 
tribution, we shall primarily focus on the case of two 
charged planar surfaces located at z = ±a with the 
charge distribution 



p (r) = aS{a 



aS(a + z). 



(3) 



We also consider a dielectric inhomogencity between the 
bounding surfaces and the ionic solution, such that the 



ionic solution is described with dielectric constant e and 
the bounding surfaces with e'. In this geometry the 
Coulomb interaction potential is composed of the direct 
interaction Uo(r, r') = (47reeo|r — r'|) _1 and the (electro- 
static) image interaction Ui m (r,r'), so that 



u(r,r') = u (r,r') 



i(r, 



(4) 



The Green's function in planar geometry can be ex- 
pressed as a sum of image charge contributions, viz. 



u(r,r') 
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.(5) 



In both sums the index n runs through negative and pos- 
itive integer values. The term corresponding to n = 
represents the direct Green's function of the particle in a 
space of a uniform dielectric constant, namely uo(r, r'). 
The unit vector k = (0,0, 1) points in the +z direction. 
The relative dielectric jump at the two bounding surfaces 
is quantified as A = (e — s')/(e + e'). For most relevant 
situations the outer permittivity is smaller than the inner 
one, e' < e, so that the dielectric jump is positive, A > 
and vanishes in the homogeneous case with £'—►£. 

The canonical partition function Zn of this Coulomb 
fluid composed of N charged particles is given by 



i J dRi • • • dRyvdu;! • • • du N x (6) 
xe^(-/?tf[R i ,h> i ;R J -,a> 3 -]) J (7) 



with the configurationally dependent ionic interaction en- 
ergy given by 

\zZ Jf drdr 'p(r; R, <*i) u(r, r') p(r'; + 
+ YlJJ drdr '/5( r ; R.i,Vi)u(r,r')p (r') + 



drdrVo(r)u(r,r')po(r'), 



(8) 



where j3 = 1/ksT. The three terms in C/[Rj, uii\ R^, utj] 
correspond to direct electrostatic interaction between 
counterions, electrostatic interactions between counteri- 
ons and fixed charges and between fixed charges on the 
walls themselves, respectively. We will furthermore make 
the standard assumption that the system is overall elec- 
troneutral, implying that the mobile countercharge ex- 
actly compensates the charge on the surfaces. 

We now proceed in the Netz [T3| fashion and perform 
the Hubbard-Stratonovitch transformation of the parti- 
tion function , where the configurational integral over 
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countcrion positions is transformed into a functional in- 
tegral over a fluctuating auxiliary electrostatic potential 
4>(r). In this way the grand-canonical partition function 
is obtained straightforwardly in the form [Io| 



z = J2 ^ N Z N 

N=0 



C / P[0(r)]e 



-0H[<t>(r)] 



(9) 



where A is the bare fugacity which is the exponential 
of the chemical potential. The prcfactor C above is the 
functional determinant of the inverse Coulomb kernel 
w _1 (r, r') while the "action" of the functional integral 
is given by 



dr 



i £ (r)e o (V0) 2 - ipo(r)^(r) 



-J du ft(r)exp \ ipjdr'p(r';r t u})(f>(r') 



.(10) 



Here we have introduced the renormalized fugacity as 
A' = A exp(|uo(r, r)), where u (r,r) is the electrostatic 
direct self-energy of a single counterions with multipolar 
charge distribution. f2(r) is the geometric characteris- 
tic function of the counterions, being equal to unity in 
the slab between the bounding surfaces and zero other- 
wise. The partition function in the above form can not 
be evaluated explicitly except in the one-dimensional case 
fl2| where the partition function evaluation is reduced to 
a solution of a Schrodinger-like equation. Nevertheless 
the field-theoretical representation of the grand canon- 
ical partition function allows one to use quite powerful 
analytical approaches that eventually lead to an explicit 
evaluation of the partition function in two well defined 
and complementary limits [l(| EH, E3 that we shall ad- 
dress later in this paper. These limits retain the relevance 
also in the case of structured counterions. 



III. DIMENSIONLESS REPRESENTATION 

One may obtain a dimensionless representation for the 
present system by rescaling all length scales with a given 
characteristic length scale. Recall that the character- 
istic distance at which two unit charges interact with 
thermal energy k B T is known as the Bjerrum length 

= eo/(47ree fcBT) (in water at room temperature, the 
value is l B ~ 0.7 nm). If the charge valency of counteri- 
ons is q then the aforementioned distance scales as q 2 i B . 
Similarly, the distance at which a counterion interacts 
with a macromolecular surface of surface charge density a 
with an energy k B T is called the Gouy-Chapman length, 
defined as \i = e$/ (2-Kq(. B o~). A competition between ion- 
ion and ion-surface interactions can thus be quantified by 
the ratio (S) of these characteristic lengths, that is 

5 = q 2 C B /n = 2irq 3 e 2 B a/e , 

which is referred to as the electrostatic coupling param- 
eter. 



In what follows, we may rescale the length scales with 
the Gouy-Chapman length, i.e., r — * r//x; hence, the 
surface separation will be rescaled as D = D/fj, or the 
rescalcd half-distance as a = a/fi. Other dimensionless 
quantities can be defined as follows. The dimensionless 
multipolar moments are defined as 



P=Po/eoq/J> and t = t Q /e qiJ> 2 
and the dimensionless pressure as 

pP 



P 



2ir£ B (a/e ) 2 



(11) 



(12) 



For simple structureless counterions the functional in- 
tegral ((9]) can be rescaled [1(3] yielding 



Z^C V[<j>(r)]e 



-H'[4>(r)]/E 



with dimensionless H'[<f>(r)]. As will be shown later an 
identical representation can be obtained also for struc- 
tured counterions, described on the level of a multipolar 
expansion of their charge density. 



IV. STRONG AND WEAK COUPLING 
DICHOTOMY 



In the absence of a general approach that would 
cover thoroughly all the regions of the parameter space 
one has to take recourse to various partial formula- 
tions that take into account only this or that facet of 
the problem [l5|. The traditional approach to these 
one-component Coulomb fluids has been the mean-field 
Poisson-Boltzmann (PB) formalism applicable at weak 
surface charges, low counterion valency and high tem- 
perature [Hldil, E3]- The limitations of this approach 
become practically important in highly-charged systems 
where counterion-mediated interactions between charged 
bodies start to deviate substantially from the mean-field 
accepted wisdom (l5j . One of the most important re- 
cent advances in this field has been the systematization 
of these non-PB effects based on the notions of weak and 
strong coupling approximations. They are based on the 
field-theoretical representation of grand canonical parti- 
tion function (flTj)) . whose behavior depends on a single di- 
mensionless coupling parameter S [H| (see below). The 
weak-coupling (WC) limit of E — > coincides exactly 
with the mean-field Poisson-Boltzmann theory [2(| and 
is based on a collective description of the counterion den- 
sity. The strong coupling (SC) limit of S — > oo is diamet- 
rically opposite and corresponds to a single particle de- 
scription. It has been pioneered by Rouzina and Bloom- 
field [HI, elaborated later by Shklovskii et al. [H, [H, [I| 
and Levin et al. [25( | , an d eventually brought into a final 
form by Netz et al. 0, GJ, GJ, ES 123, S3 • 

These two limits are distinguished by the pertaining 
values of the coupling parameter S which can be intro- 
duced in the following way. 
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The regime of S <C 1 is the case of low countcrionic va- 
lency and/or weakly charged surfaces, and is referred to 
as the weak-coupling limit. It is characterized by the fact 
that the width of the counterion layer p is much larger 
than the separation between two neighbouring counte- 
rions in solution and thus the counterion layer behaves 
basically as a three-dimensional gas. Each counterion 
in this case interacts with many others and the collec- 
tive mean-field approach of the Poisson-Boltzmann type 
is completely justified. 

On the other hand in the strong coupling regime S ^> 
1, which is true for high valency of counterions and/or 
highly charged surfaces. In t his cas e the mean distance 
between counterions, a_i_ ~ yf eoq/a, is much larger than 
the layer width (i.e., a±/fi ~ 1), indicating that 

the counterions arc highly localized laterally and form 
a strongly correlated quasi-two-dimensional layer next 
to a charged surface. In this case, the weak-coupling 
approach breaks down due to strong counterion-surfacc 
and counterion-counterion correlations. Since counteri- 
ons can move almost independently from the others along 
the direction perpendicular to the surface, the collective 
many-body effects that enable a mean-field description 
are absent, necessitating a complementary SC descrip- 
tion EiEiBl. 

These two approximations allow for an explicit and ex- 
act treatment of charged systems at two disjoint limiting 
conditions whereas the parameter space in between can 
be analysed only approximately and is mostly accessible 
solely via computer simulations. As will become clear 
in what follows the WC and SC limits remain valid also 
for structured counterions described with a multipolar 
charge distribution. 



V. MEAN-FIELD LIMIT 

The mean- field limit, 5 — > 0, is defined via the saddle- 
point configuration of the Hamiltonian IjlOp [20| and is 
valid for a weakly charged system. This leads to the fol- 
lowing generalization of the Poisson-Boltzmann equation 

££ VVo(r) = -A'// dr' du p(r; r', w) fi(r') X 



x exp 



-l3jdr"p(r";r',u)Mr") 



(13) 



for the real- valued potential field ipo — —i(f>o. Note that 
the dielectric discontinuity at the boundaries is irrele- 
vant within the mean-field theory [2^ | since in the planar 
geometry this is effectively one-dimensional theory. Tak- 
ing now the counterion density function as a sum of the 
monopolar, dipolar and quadrupolar terms, where the 
mean- field depends only on the transverse coordinate z, 
the above equation can be written in dimcnsionlcss form 



ip" = ~ I dx 0(x, z) (u(z) - pxu'(z) + tx 2 u"(z)) 



<x pi{z) + p 2 {z) + p 3 {z), 



(14) 



where we have defined dimcnsionlcss potential ip = 
fleoqipo with corresponding derivatives ip' = (3eoqpiJj' and 
tp" = fteoqp ipQ and dimcnsionlcss multipolar moments 
p and t. We defined the orientational variable x = cos 9, 
where 9 is the angle between the z-axis, assume to coin- 
cide with the normal to the bounding surfaces, the direc- 
tor of the uni-axial counterion is n and the integral over 
this variable gives the orientational average. We have 
also introduced 

u(z) = C exp (— tj) — pxtfj' — tx 2 ip") , 

which is obviously the local orientationally dependent 
number density of the counterions and p(z) is the corre- 
sponding orientationally averaged charge density of the 
counterions. In the case of higher multipoles the num- 
ber density and the charge density of the counterions 
are not proportional. Expressions pi(z) are simply the 
orientationally averaged multipolar charge densities, i.e., 
i = l for monopolar charge, i = 2 for dipolar charge etc., 
that are simply proportional to the three terms in the 
integrand of (|14[) . The corresponding prefactor is simply 
obtained by appropriate normalisation of pi to satisfy 
electro-neutrality condition of the system. 

The above PB equation has to be supplemented by an 
appropriate boundary condition at z = ±a corresponding 
to the electroncutrality of the system by taking into ac- 
count surface charges, a. The constant C is set by these 
boundary conditions, as is the case in the standard PB 
theory. In dimcnsionless units these boundary conditions 
read ip'(±a) = =f2. Obviously in the case of counterions 
with only monopolar charge distribution the above set 
of equations reduces to the standard Poisson-Boltzmann 
theory. 

The characteristic function &l(x, z) for the parallel 
plane geometry simply excludes the counterion config- 
urations that would penetrate the bounding walls and 
depends on the geometric form of the counterions. Since 
in what follows we will not be interested in steric effects, 
we will assume that all counterions are pointlike, and so 
Q,{x,z) = 1. Steric effects have been studied elsewhere 



Because of the similarity with the standard PB equa- 
tion for monopolar charges one is led to believe that the 
pressure in an inhomogeneous system of multipolar coun- 
tercharge can be derived in the same way as in the stan- 
dard PB theory, via the s.c. contact value theorem f30j ] . 
Indeed this can be proven exactly, namely the mean-field 
pressure can be obtained from the first integral of the PB 
equation (|14[) . After some manipulations its first integral 
is obtained in the form 

P= -\ip' 2 + \ J Ax [u + pxutp' + tx 2 (wp" - • 

(15) 
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Here P is the dimensionlcss equilibrium pressure in the 
system, defined as ([T2"jl . The r.h.s. of the above identity 
can be calculated at any arbitrary point \zq\ < a as its 
value is independent of zq. The choice of the actual point 
is governed by the symmetry of the system, as is usual 
also in the standard PB theory. Since we delve only on 
the symmetrical solution, with inversion symmetry cen- 
tered on z — 0, the derivative of the mean potential at 
the mid-point must vanish, i.e., ip'{zo = 0) = 0. Thus in 
this case 



p = 0, which then reads 



P 



dxu(0)(l-tx 2 \i/j"(0)\) 



(16) 



While the first - van't Hoff - term corresponds to repulsive 
interactions of the standard PB type, the second term in 
this particular geometry and symmetry entails attractive 
interactions between the bounding surfaces. Note that 
V>" < 0. 

Let us take a closer look at the above expression. For a 
polyelectrolyte chain, i.e., an extended and flexible coun- 
terion, one can derive a similar type of pressure formula 
[3lj . with a negative term that contributes an attractive 
part to the force equilibrium. This attractive part in the 
case of polyelectrolytes is due to polyelectrolyte bridg- 
ing interactions [HJ that stem from the connectivity and 
flexibility of the polyelectrolyte chain. 

Does equation (| 16[) have a similar physical content? 
This interpretation certainly does not seem likely for 
point countcrions, for which the above pressure formula 
was derived. The structure of the first integral of the PB 
equation (|16|) seems to be saying that apart from the ideal 
contribution to the equilibrium pressure, a term propor- 
tional to the midpoint number density, one also finds 
an electrostatic contribution that is due to the interac- 
tion between the monopolar and the quadrupolar part of 
the counterion charge density across the midplanc and is 
proportional to the square of the midpoint density, see 
(fT4]l . The attractive part thus does not look like bridging 
which has its origin in the connectivity and flexibility of 
the polyelectrolyte chain, but more like a virial expansion 
in terms of the multipolar interactions. This is of course 
only true for point-like counterions. For extended coun- 
tcrions a bridging interpretation would be more appropri- 
ate as was already clear in the early studies of structured 
counterions 13411. 



A. First order in t 

The above PB equation is a fourth order highly non- 
linear intcgro-diffcrential equation and therefore difficult 
to handle with conventional numerical procedures. Since 
here we are not particularly interested in the weak cou- 
pling results, we want only to show that the effect of 
quadrupolar moments, t, on the interaction pressure is 
to induce a small attractive contribution. In order to 
show this we expand (|14[) to the first order in t and put 



V/' = -Ce"*(l + \t^' 2 - §#"). 



(17) 



This equation can be solved much more easily with con- 
ventional numerical methods. Although we must be 
aware that it is valid only for small t. The correspond- 
ing formula for interaction pressure in the 1st order in t 
expansion reads 



P = iCe" 



Ko) 



(18) 



Due to a convenient cancellation the formula contains no 
explicit t-dependence, but the pressure still depends on 
t implicitly, due to the t dependent potential ip. As can 
be shown numerically for small t the correction to the 
standard PB pressure depends linearly on t. 



B. Numerical results in the mean-field limit 

Here we are not interested in the details of the mean- 
field results-they were analyzed in detail before [1, 0, [1]- 
but list them for completeness anyhow. In what follows 
we delimit ourselves to counterions that posses monopo- 
lar and quadrupolar charge, so we do not take into ac- 
count any dipolcs (p = 0). 

The density profile on Fig. [2ja) corresponds to 
monopolar density p\ (z) , (|14p . Counterions are localized 
preferentially in the vicinity of both surfaces due to mu- 
tual repulsion that prevents them from being localized at 
the center. What we then observe is that with increase 
of the quadrupolar moment, t, countcrions concentrate 
at surfaces even more. This can be easily explained. The 
potential energy of every quadrupolar particle in electro- 
static potential is 



PW = tip" cos 2 



(19) 



Since the second derivative of the mean-field potential ip" 
in the symmetric case considered here is a concave func- 
tion of the coordinate z, this means that the quadrupolar 
force, F = —dW/dz, acts away from the center toward 
both surfaces. 

The insight into the orientation of quadrupolar mo- 
ments can be obtained through the average of the square 
of variable x, namely (x 2 ), 



(x 2 ) 



/ x 2 u(x)dx 



J u(x)da 



(20) 



It is instructive to define an orientational order parameter 

S=^{3(x 2 )-l). (21) 

This order parameter can have values in the range 
from —1/2 to 1. In a totally disordered system where 
quadrupoles point randomly in all possible directions 
(x 2 ^ = 1/3, and the orientational order parameter van- 
ishes, S = 0. In the completely ordered case, 5 = 1 and 
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FIG. 2: Mean-field results for various t in first-order approximation, (|17|l , (a) Counterion density profile, 1)14) 1 . By increasing 
i in the mean-field limit counterions are depleted from the center toward both surfaces. All the densities are rescaled to 
correspond to the same area under the curve, (b) Order parameter profile, evaluated from (|20)) and 1)21) 1. (c) Dimensionless 
pressure-distance curves obtained from the first-order approximation of (fl7)l and l[I5)l. Quadrupolar contributions are attractive. 



(x 2 ^ = 1, so that all quadrupolar moments are paral- 
lel with the z-axis and perpendicular to the walls, 6 = 0. 
The other extremum, S = —1/2, corresponds to (a; 2 ) = 
so that all quadrupolar moments are perpendicular to z- 
axis and parallel with the walls, 6 = w/2. 

Fig-IHb) shows the numerically obtained profile of the 
orientational order parameter. It increases with increas- 
ing quadrupolar strength t and is larger at both surfaces 
than at the center. The order parameter S > indicates 
that quadrupolar moments are preferentially aligned par- 
allel to the z-axis. According to (| 19[) the electrostatic en- 
ergy can be minimized by increasing cos 2 9, since tp" < 0. 
So (cos 2 6*) really is preferably increased above the value 
of 1/3. 

The quadrupolar contribution to the pressure for t > 
turns out to be to be attractive, Fig. ^ as was shown 
before @, 0, HI- Again this attractive contribution is 
obtained here on the mean-field level where the inter- 
counterionic correlations, essential in the strong coupling 
limit, play no role. It is indeed the intra- counterion corre- 
lations (entering here via the rigid structure assumed for 
counterions) that lead to such attractive contributions. 
At this juncture it does not make much sense to us in 
continuing the multipolar expansion to yet higher orders, 
that would become relevant for even more aspherical and 
elongated charge distributions that would obviously en- 
tail also some molecular flexibility. In that case one could 
use the well developed theory of polyelectrolyte mediated 
interactions with much better confidence 13511. 



VI. STRONG COUPLING LIMIT 

In the strong coupling limit, 5 > 1, the system is 
highly charged and counterions are highly correlated. 
The partition function §§§ can be approximated by the 
two lowest order terms in the virial exp ansion with re- 
spect to the (renormalizcd) fugacity [l0( 



AO) 



0(X' 2 ). 



The second term corresponds to a one-particle partition 
function as the SC limit is effectively a single particle 
theory. Z^l and are then given by 



^SC = CX P 



1/3/ drdrV (r) U (r,r'Vo(r') 



(23) 



which is the exponential of the interaction between bare 
surface charges, and 

7 {l) 



y // drdr'p(r;R,w)u(r,r')/5(r';R,u;) 



(24) 



which is the partition function of all possible (single) 
counterion configurations. The interaction potential in 
this part of the partition function is composed of the di- 
rect and image electrostatic interactions (U]). 

In the case of two charged surfaces with uniformly 
smeared surface charge density ((31), the surface charge 
electrostatic potential does not depend on the z- 
coordinate; it is spatially homogeneous and is given by 



/ drXr,r')po(rO = - — . 

J ££o 



(25) 



The corresponding potential energy of a counterion in 
this surface electrostatic potential is given by 



J J drdr'p(r; R, u))u(r, r')p (r') 



ee 



-e q. 



(26) 



(22) 



Since all the terms in density operator (fT]), except the 
first one, depend on the gradients, i.e., spatial deriva- 
tives, the counterion energy in a homogeneous external 
electrostatic potential depends only on the first, monopo- 
lar term (f2"6")) . That means that higher multipoles do not 
interact directly with planar surface charge, but we must 
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emphasize that this is not the case in inhomogeneous po- 
tential at curved surfaces, e.g., at cylindrical or spherical 
surfaces or inhomogencously charged surfaces. 

As for the self-energy contribution, the second term in 
the exponent of (|24|) , it only picks up terms from the z- 
dependent parts of the image self- interaction, u; m . The 
direct self-interaction, ito, docs not depend on coordi- 
nates so it can be discarded. The self-image energy is 
then 



action between a monopole and its own monopole im- 
age, therefore it is proportional to q 2 . The second term, 
1U2, corresponds to the monopolc-dipole image interac- 
tion which is obviously the same as the dipolc-monopolc 
image interaction, w^, and so on for all the higher order 
terms. 

These image self-interaction terms can be expressed in 
a dimcnsionlcss form as 



drdr'/5(r; R, w)u im (r, r')p(r'; R, lj) — 'S^vjj{z, cosfl). 

i=l 

(27) 

Self-image contributions are among monopolar, dipolar 
and quadrupolar moments of the countcrions interacting 
with their own electrostatic images of monopolar, dipolar 
and quadrupolar moments. Summing up all these con- 
tributions we remain with nine terms of which only six 
are different. Writing them up in extenso we obtain 



Wi(z,cos8) — I; (3 Wi(z, cos 9). 



(29) 



wi(z 


cos 9 


) = q 2 el u im (R,R), 


w 2 {z 


cos 9 


) = qe p (n ■ V)u im (r, R) r=R , 


w 3 (z 


cos 9 


) = qe t (n ■ V) 2 w im (r, R) r=R , 


w 4 (z 


cos 9 


) = w 2 (z, cos#), 


w 5 (z 


cos 9 


)=PI (n- V)(n- V')uim(r,r')| r , =r=R , 


w 6 {z 


cos 9 


) = Po t Q (n • V)(n • V') 2 Wi m (r,r')| r/=r=R 


w 7 {z 


cos 9 


) = wa(z, cos 9), 


w 8 (z 


cos 9 


) = wq(z, cos 9), 


Wg(z 


cos 9 


) = t 2 (n.V) 2 (n.V') 2 Uim (r,r')| r , =r=R . 



(28) 

Here, the first expression, W\, corresponds to the inter- 



The above expressions were derived for the planar geom- 
etry. Similar expressions but with different image poten- 
tials can be derived also in other geometries with dielec- 
tric inhomogencitics. Note that if there is no dielectric 
discontinuity, so that u; m (r,r') = 0, they are all identi- 
cally zero, Wi=0\ 

Thus one can conclude at this point that in the dielec- 
trically homogeneous case the higher order multipoles are 
completely irrelevant on the SC level. Only the monopo- 
lar term survives in the partition function. This is prob- 
ably one of the most important conclusions of this work, 
so let us reiterate it: in a dielectrically homogeneous case 
in plan-parallel geometry the structure of the countcrions 
as codified by their multipolar moments plays absolutely 
no role in the strong coupling limit! 

For our case of two charged planar surfaces we eval- 
uate the above expressions by using the Green's func- 
tion ([5]); which was derived for planar geometry with 
a step-function dielectric profile at the two boundaries 
z = ±a. We get the following somewhat cumbersome 
expressions for the dimcnsionlcss self-image interactions 
defined above 



J 



"'2 



Wi[Z,Xj = — £1 > , ,, 

K ' J 2 ^ (no) 2 - 

n odd v ' 



{z,x) = ^-Epxz 



naA n 



odd 



[(na) 2 — z 2 ] 2 ' 



1 A 

^ n even 



E 



. (na) 2 + 35 2 . _ 

™ - ~-2 _ =213 



(na) 3 ^— ' \(naY — z 

K ' n odd LV ; 



~ ,~ ^ U 2 f/, o 2x A " 2>V- - (n&) 2 +3S 2 

w 5 (z, x) = -^p 2 I (1 - 3x 2 ) }2 7^3 +0- + x ) E na ±^ — A 

^ n ever 



n odd 



we(z, x) — —Spt(l + x 2 )x z 
w 9 (S,x) 



(na) 

_ (na) 2 + z 2 
naTT — 6i ^rrA 11 



[(t^) 2 - z 2 } 3 



[(nhf 



32 



n odd 



[(nay 



(30) 



3 , \ , , v-^ A " 9 dl (na) 4 + I0(na5) 2 + 5z 4 
Et 2 { (3 - 30x 2 + 35x 4 ) 7-^5 + ( 3 + 2x + 3x ) E ™ — rT — Jz^r- A" 



All indices n here run only through positive values. Note again that self-image contributions, u),-, vanish when 
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A = 0. Taking all this into account, the strong cou- 
pling interaction free energy, (3F = — k-QTlnZ^,, i.e.. 
the part of the free energy that depends on the inter- 
surfacc separation, can be written in dimensionless form 
as 



F/i = 25-21n 



ix expl 



i=l 



(31) 

Just as in the WC case we have again assumed that 
the counterions are point-like particles so that the char- 
acteristic function f2 does not depend on their coordinate 
z and therefore the integration goes from —a to a. This 
means that we also disregard the possible entropic ef- 
fects due to the finite size and anisotropy in the shape of 
counterions. These entropic contributions to the parti- 
tion function are relevant only for intersurface separation 
on the order of the size of the countcrion or smaller. In 
that regime of separations other, much stronger effects 
would come into play and compete with ionic finite size 
effects, thus these type of effects are not the focus of this 
paper. 

The corresponding dimensionless pressure is simply ob- 
tained by taking the derivative of the free energy (j3"T]) 
with respect to wall separation, D = 2a, 



P = - 



dF/A 
2 da ' 



(32) 



It is again obvious from here that in the dielectrically 
homogeneous case, i.e., the case with no electrostatic im- 
ages, where Wi(z, x) = 0, the strong coupling limit is 
given exactly by the monopolar term (the first term in 
(|24p). Thus without the images we remain with the same 
form of the strong coupling interaction free energy as in 
the case of monopolar point charges. It is given by 



F/A = 2a - 21n< 
with the corresponding pressure 



(33) 



P 



1. 



The higher order multipoles thus make no direct con- 
tribution to the strong coupling interaction free energy 
or forces between the bounding charged surfaces without 
dielectric inhomogencitics. 

This is a very powerful result whose ramifications in 
fact impose rather stringent limits on the significance of 
the multipolar expansion of counterionic charge. It ap- 
pears that for highly charged counterions, that are in 
fact the only ones where the multipolar expansion really 
makes sense, most of the electrostatics is properly cap- 
tured by the monopolar term. Higher order multipoles 
simply do not contribute to the pressure in the system 
in the SC limit. Of course all of this is valid in the limit 
of homogeneous dielectric properties for planar surfaces 
without any image effects. 



The countcrion density profile can be extracted from 
(|3ip as an integrand in the second term, 



,1 9 

p(z) oc / dx cxp( — Wj(z, : 
i=i 



(34) 



In the case without the dielectric mismatch the density 
profile is simply homogeneous, p(z) = const. Dielectric 
images induce an additional repulsion between the coun- 
terions and the surface charges pushing them towards 
the midplane region between the two charged bounding 
surfaces. 

Oricntational order of quadrupolcs can be inspected 
via averaged (x 2 ) which is here defined as 

f x 2 exp(-^Li Wi(z,x))dx 

(*-) = - . ./ ■ (35) 



/ exp(^- X\ =1 Wi(z, x)\ Ax 

We can now use the same definition of order parameter 
S as in WC, (HQ). 



A. Numerical results in the SC limit 

We now present some numerical results for the SC 
limit. Without any dielectric mismatches the density pro- 
file of the monopolar counterions is homogeneous, i. e., it 
does not depend on z, which is a well known result [10| . 
In the case of multipoles with dielectric images this is no 
longer true, but the counterion density does remain an 
even function of z. 

On Fig. [3] we show the numerical results for coupling 
parameter strength S = 50, dielectric jump A = 0.95 
and 5 = 1. Contrary to the WC situation, here counte- 
rions are localized in the vicinity of the midpoint be- 
tween the two bounding surfaces due to image repul- 
sion. Even without quadrupolar contributions (t = 0) the 
monopolar interaction, w\, is repulsive and long ranged, 
u>i ~ z~ x [2S| . compressing the counterions to the mid- 
point. 

When the quadrupolar parameter t, is increased a bit 
the midpoint peak first starts to widen. This is due to the 
attractive monopole-quadrupole image contribution, rep- 
resented by the W3 and W7 terms which are proportional 
to t. As t increases further (approximately for t > 8a 2 /45 
for large AS) the repulsive quadrupolc-quadrupolc im- 
age contribution, wg, which is proportional with t 2 , takes 
over and confines counterions to the center even more so 
that the density peak sharpens up. 

For the value AS large enough, so that density pro- 
file becomes bell-shaped, the SC density profile can be 
approximated by a gaussian curve with a variance, i.e., 
peak width, of the form 



Sz 2 



1 



AS 5 4 



3a 2 t + ±f i 2 ' 



(36) 
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(c) 

rious t. (a) Counterion density profile, (b) Order parameter 
profile. Counterions are less ordered at the center. Order is increased with increasing t. (c) Quadrupolar pressure contribution, 
AP = P t — Pt=o- Contribution is repulsive at small separations and becomes attractive at larger separations. At very large 
separations it goes to zero. 



As can be seen, small t widens the density peak, but fur- 
ther increasing t over the value t > 8a 2 /45 the peak nar- 
rows. Also by increasing parameters S, A the midpoint 
density peak narrows. This also means that the repul- 
sive pressure between the walls increases. As the density 
peak narrows and the density of the counterions at the 
midpoint plane increases, the counterions become more 
ordered in the quasi two-dimensional midpoint layer, im- 
plying also that the correlations are stronger and there- 
fore the SC limit should a fortiori be even more appro- 
priate. 

Fig. [3^b) shows also the orientational order parameter 
for the SC limit. Contrary to the WC results, here the 
orientational order parameter is negative, implying that 
quadrupoles arc preferably aligned perpendicular to z- 
axis. The parameter reaches its extremum —1/2 at both 
walls where quadrupoles are perfectly perpendicular to 
z-axis and parallel to the walls. This is due to the strong 
quadrupole-quadrupole image repulsion, wg, which has 
its minimum at x — and is not a purel y st eric effect as 
in simulations of long charged stiff rods [36| • 

Since monopolc-monopole image interaction, wi, has 
the longest range, we can expect that for large separa- 
tions between the walls, 25,, this interaction will dominate 
and higher multipoles can be again neglected. Fig. (3{c) 
shows the pressure difference AP, which is the differ- 
ence between the pressure with finite t and the pres- 
sure without the quadrupolar contributions, t = 0. As 
expected this difference indeed goes to zero at large 
separations. At moderate separations the quadrupo- 
lar contribution becomes negative, which is caused by 
monopolc-quadrupolc image attractions, W3 and W7, that 
are shorter ranged than wi. At even smaller separa- 
tions quadrupole-quadrupole image repulsion, wg, be- 
comes dominant, since it has the shortest range and 
varies approximately as ~ a~ 5 . The larger the t the 
smaller is thus the maximal quadrupolar attraction at 
larger distance since it is overwhelmed by repulsive con- 
tributions. 



VII. MULTIPOLES, CORRELATIONS, 
BRIDGING AND ALL THAT JAZZ 

There appear some similarities between the polyelec- 
trolyte bridging interaction (32l . |33| and the attraction 
seen with strongly charged counterions. Recent simula- 
tions and density functional results (36j in fact do accen- 
tuate a close connection between polyelectrolyte bridg- 
ing and strong-coupling electrostatics, a line of reasoning 
that we will explore here in finer detail. 

Turesson et al. indeed find out that for flexible 
polyolcct.rolyt.os bridging attraction appears to be the 
dominant source of attractive interaction. On the other 
hand, for stiff charged rods the strong adsorption to the 
bounding charged walls prevents bridges to form. In- 
stead a very strong correlation attraction is apparent at 
shorter separations. Flexibility thus extends the range 
of attractive interactions as their nature changes from 
short range correlation attraction to longer ranged poly- 
electrolyte bridging. Similar conclusions have been also 
reached in [37j . 

Here wc will formulate the exact correspondence be- 
tween the intuitive notion of "bridging" in the case of 
multivalent counterions and the definition of "bridging 
interaction" in the case of charged polymers. Let us 
analyse first the partition function for two disparate 
systems: a counterion Coulomb fluid confined between 
two strongly charged walls, and a flexible polyelectrolyte 
chain confined in the same geometry with weakly charged 
walls. For point counterions between two charged sur- 
faces the grand canonical partition function in the SC 
limit [38| can be written to the first order in fugacity as 

Z = zjV+\' £&> + ..., (37) 

where 

4°C=C I V[<Kt)] exp(-i/fe J £(r)(V<Mr)) 2 dr + 

+ i(3 / Po (r)0(r)drY (38) 
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and 



Z$ =C J dR Q(R)Jvy>{r)] exp (- ±/3e Je(r) ( V(^(r)) 2 dr + 
+ i/3 J p o (r)0(r)dr + i(3qe J p(v - R)0(r)drj . (39) 



The above partition function is of course exactly the same 
as the one in (|2"2"|) after one evaluates the functional in- 
tegrals indicated above explicitly. Here po(r) is again 
the external fixed charge on the two bounding surfaces 
located at z = ±# and the shorthand p(r — R) is intro- 
duced for p(r - R) = 5 3 (r - R). 

On the other hand the partition function for a flexible 
polyelectrolyte chain is given by [39l | 

Z P =cJv[0(r)} exp^-i / 3£ o y £ (r)(V0(r)) 2 dr+ (40) 

+ i/3 Jp o (r)0(r)dr + In JJ ^(r, r'; TV) drdr'^j . 

where C^(r,r') is the Green function of a polyelectrolyte 
chain in external electrostatic field </> given by [40| 



<r'=R(iV) 

g 4> {v,r'-N)= V[R(n)]x 

Jr=R(0) 

3 f N fdR(ny- 



x exp I 



y dri 



dn + i/3geo 



(41) 

jf 0(R(n))dnY 



where R(n) is the coordinate of a polyelectrolyte 
segment n and the boundary condition is set as 
limjv-»o G<t>{ v i r '; N) = S 3 (r — r'). Representing the 
Green function via a corresponding Edwards equa- 
tion one can derive that the polyelectrolyte monomer 
density p^v^N) is given by a functional derivative 
of In J" G$(r, r'; N) drdr' with respect to the fluctuat- 
ing electrostatic potential </>(r) [4(|. The polyelec- 
trolyte monomer density has to satisfy the condition 
limAr^o/O0=o(r;iV) = <5 3 (r). 

In the planar geometry with two bounding surfaces, 
one can center the monomer density at the origin of the 
coordinate system, whose z axis is directed along the 
surface normal and has the origin at the midplanc. De- 
noting R = (0,0,0) one can thus write p</,=o(r; AT) = 
P0=o(r — Ro; N). In the limit of small electrostatic po- 
tentials, the lowest order (i.e., the first-order) solution of 
the polyelectrolyte Green function can be written as 



In JJ g^r, r';N) drdr' 



i(3qe Q J ^=o(r-Ro; N)<f>(r)dr. 

(42) 

By construction this solution clearly corresponds to a 
weak coupling of the polyelectrolyte to the external field 
generated by the fixed charges at the boundary of the 
system. For large values of this charge higher orders in 
the expansion (|4"2")l would certainly have to be taken into 
account. 



In the case of weak coupling the partition function of 
a polyelectrolyte chain is then given approximately by 



Zp~4^ =o (r-R o ;A0] = 



(43) 



-C J Vy>(r)] cxp^-i/feo J e(r)(V^(r)) 2 dr + 



i/3 / p (r)<Kr)dr + i(3qe / p 0=o (r - R ; N)</>(r)dr 



If we now compare the expressions (|39[) and (|43[) we can 
establish the following identity 



Znl! = lim 



dRn(R)Z^j^ =0 (r-R;JV)], (44) 



with the simple (first order SC) counterion partition func- 
tion on the l.h.s. and the polyelectrolyte (zeroth order 
WC) partition function on the r.h.s. This "duality" re- 
lation connects the partition function of a strongly cou- 
pled point counterion system with the partition function 
of a weakly coupled flexible polyelectrolyte system. Since 
Z$c f° r finite N describes polyelectrolyte bridging inter- 
actions it is clear that the above formula establishes a for- 
mal connection between strong coupling point counterion 
attractive interactions and weak coupling polyelectrolyte 
bridging interactions. 



VIII. CONCLUSIONS 

Using field-theoretic methods we derived a description 
of the counterion- mediated electrostatic interaction when 
the countcrions posses internal degrees of freedom such as 
a rotational axis. We concentrate on the symmetrically 
charged planar surfaces with a dielectric mismatch on 
both sides and countcrions with quadrupolar moments in 
between. We analyse this system in two different regimes, 
namely in the weak coupling (mean-field) and the strong 
coupling limit. The former one describes the case of low 
surface charge and low counterion valency, whereas the 
latter one describes high surface charge and high counte- 
rion valency. 

In the mean-field limit, we derived the Poisson- 
Boltzmann equation for counterion density in the case 
that the counterions posses dipolar as well as quadrupo- 
lar moments. Since the equation is highly non-linear, 
we solved it numerically only in the first order approx- 
imation, which should be valid for small quadrupolar 
moments, t. As is already known the dielectric mis- 
matches play no role in the mean-field theory and have 
therefore no effect on quadrupoles [29| . One could add 
those corrections in by hand [4l|, but such an approach 
does not strictly corresponds to the mean- field analysis. 
Quadrupolar interaction affects the counterion density 
distribution as well as interactions between the charged 
surfaces. Counterions are depleted from the central re- 
gion and concentrated in the vicinity of both surfaces. 
The orientation of quadrupoles is preferentially parallel 
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to the z-axis as the moment t is increased. The highest 
alignment is reached right next to both bounding sur- 
faces. The quadrupolar contribution to the interaction 
pressure is attractive. 

In the strong coupling limit higher order multipoles 
also play an important role in the intersurface interac- 
tions but only for a diclcctrically inhomogcneous case. 
Without any dielectric discontinuities in the system the 
higher multipoles simply do not matter in plan-parallel 
case and the monopolar part of the partition function 
captures all the strong coupling effects. This is a com- 
pletely general result for plan-parallel surfaces and should 
be of some importance in assessing the electrostatically 
mediated forces between strongly charged planar sur- 
faces. As already mentioned the multipoles can have 
effects near curved surfaces even without dielectric dis- 
continuities. If compared to the case of monopolar coun- 
terions, higher multipoles interact only via the orien- 
tationally dependent part of the image self-interaction. 
In our analysis we focused on the counterions with 
monopolar and quadrupolar moments, so that the in- 
teraction is composed of three contributions. These are 
the monopole-monopolc image interaction, which is long 
ranged and repulsive, the monopolc-quadrupole image as 
well as quadrupole-monopolc image interactions that are 
shorter range and attractive. Finally, very short ranged 
quadrupole-quadrupole image contribution is repulsive 
and plays a dominant role at very small distances. 

We are at present involved in extensive Monte Carlo 



simulations in order to explore the validity of the for- 
mal developments described above. Preliminary results 
completely vindicate the theoretical analysis. 
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